Morphological regions and oblique incidence dot formation in a model of surface 

sputtering 
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We study solid surface morphology created by off-normal ion-beam sputtering with an atomistic, 
solid-on-solid model of sputter erosion. With respect to an earlier version of the model, we extend 
this model with the inclusion of lateral erosion. Using the 2-dimensional structure factor, we found 
an upper bound fj, ~ 2, in the lateral straggle /x, for clear ripple formation. Above this upper 
bound, for longitudinal straggle a > 1.7, we found the possibility of dot formation (without sample 
rotation). Moreover, a temporal crossover from a hole topography to ripple topography with the 
same value of collision cascade parameters was found. Finally, a scaling analysis of the roughness, 
using the consecutive gradient approach, yields the growth exponents /3 = 0.33 and 0.67 for two 
different topographic regimes. 

PACS numbers: 05.10.-a,68.35.-p,79.20.-m 



I. INTRODUCTION 

Low energy bombardment of the surface of amorphous 
materials Q, semiconductors (amorphized by the sput- 
tering processs) 0- S C3] ■ and metallic materials (at low 
temperature) |5j, by a beam of ions at off- normal inci- 
dence, often lead to ripple pattern formation. The rip- 
ple orientation is perpendicular to the projection of the 
ion beam direction, onto the surface plane, for small inci- 
dence angles, and parallel to the projection for grazing in- 
cidence. However, for a metallic surface with anisotropic 
diffusion, the orientation is perpendicular to a crystallo- 
graphic direction (the one favored for diffusion) at small 
incidence angle Q- 

According to the continuum theory 0,0, ripples arise, 
for all angles of incidence 9, from the curvature depen- 
dence of the sputter yield. The basis of the calculation 
is the Sigmund distribution describing the energy 
deposited by the incoming ion. This distribution results 
from a study of the collision cascades created by the pene- 
trating ion. The Sigmund distribution is parametrized by 
the depth a, the longitudinal straggle <r, and lateral strag- 
gle fj, of the energy, cf. 0J. The wavelength A obtained 
within the continuum theory is given by A = 2wy/2K/\v\, 
where K is the surface diffusivity, and v is a negative sur- 
face tension coefficient, the latter depending on a,/z, a 
and 9. Typical wavelengths are of the order of tens of 
nanometers. In the absence of non-linearities, the 9 de- 
pendent v, being different along the parallel and per- 
pendicular directions to ion projection, govern the ripple 
orientation. However, contrary to the predictions of the 
continuum theory Jccf. Figs. 12 and 15 of ||), no ripples 
were observed in Q for 9 < 40° under Xe ion irradiation 
of Si. 

Other studies have shown that this influence of col- 
lision cascading on topography evolution extends to 
quantum dot formation at 9 = on non-rotated 
substrates^, and at 9 > on rotated substrates [Tlj . 
Also, collision cascade statistics has recently been found 



to shift \12\ the Sigmund energy distribution away from 
the Gaussian which has been the basis of the theoretical 
models so far in use [71, [lj, [l5| . But till now the ripple 
phase boundaries arising from such influence are thought 
to be time (i.e. ffuence) independent. Furthermore, up 
till now, dependence of ripple formation on the lateral 
straggle //, though acknowledged, is often ignored on ac- 
counts of isotropic or symmetric cascading for simplicity. 

Recently a discrete MC model of sputter-erosion was 
introduced [l^ which, like the continuum theories, is 
based on Sigmund's spatial distribution of the kinetic 
energy transfered by an impinging ion. Its very recent 
application to the study of ripple motion |16| indicates 
its capture of the universal features of material-surface 
modification by ion beam treatment, and affords us a 
way of exploring the different phases in the surface mor- 
phology. 

In this paper, we extend this model with inclusion of 
lateral erosion. This includes stronger non-linearities in 
the sputtering process. Hence, the exponential-growth 
tendency of the ripple amplitude occuring within the lin- 
ear theory is more stabilized [ItI llS|. Furthermore, we 
show that the angle 9 r above which ripple formation oc- 
curs depend on the longitudinal straggle a of the ion 
beam, with higher a resulting in lower 9 r . Moreover, 
we find that there are two distinct creation mechanism 
for ripple: In one region of the parameter space, there 
is a transition from hole topography at an early time, 
to ripple topography at a later time. For other param- 
eters, the ripples are created from a simple rough to- 
pography at early times. Our results also indicate the 
crucial role of the value of (i for ripple formation. If we 
chose n larger enough, then we find dots/nano-sized is- 
lands similar to those observed by 0, |2(j , and predicted 
by 0, EI without sample rotation. Our results sug- 
gest that by using new projectil/target combinations one 
might find yet unobserved surface topologies. We men- 
tion below some promising projectil/target combinations 
which might serve as guidelines for experimentalists. 
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In the next section, we describe our simulation model. 
In the main part we present and discuss our results. We 
finish with a summary and discussion. 

II. NUMERICAL SIMULATION METHODS 




FIG. 1: Sketch showing lateral erosion of surface atoms, and 
collision cascade geometry, as described in the text. Erosion 
of atom at position I or 2 triggers a relaxation, which ensures 
there is no overhang; the atom at position 3 or 4 is only 
eroded, without surface relaxation. 

We study system sizes L x L with periodic boundary 
conditions. The material surface is defined by a time- 
dependent discrete height function h(x,y,t) [(2 + 1)D 
solid-on-solid (SOS) model] which is initially flat, i.c, 
h(x,y,0) =constant. Hence, overhangs cannot occur. We 
have included surface roughening via sputtering and sur- 
face smoothing via surface diffusion. Each simulation 
step consist of the sputtering generated by one ion and a 
certain number of diffusion steps. The simulation time is 
measured in terms of these steps and corresponds to the 
fluence in experiments. For details see Ref. [jij- Here we 
only display the main ingredients of the methods, for the 
convenience of the reader, and describe the extensions 
with respect to the original method. 

We simulate the sputtering process as a combination 
of ion motion and erosion of atoms. The ion source, at 
some random position in a plane directly above the sur- 
face plane, projects the ions along a trajectory inclined at 
9 to the vertical and with azimuth </>. After penetrating 
the solid through a depth a, measured along the ion tra- 
jectory, each ion comes to rest and distributes its kinetic 
energy e. The erosion process is governed by the gener- 
ally accepted Sigmund fl3l ] distribution. A surface atom 
is eroded with a probability proportional to the energy 
deposited there, which is given by Sigmunds's formula 

e ( d l d 2 , 1 



tory. dj_ is the distance of the surface atom perpendicular 
to the ion trajectory. Hence, we use the local coordinate 
system of the ion with origin at the point of penetration, 
and the z-axis coinciding with the ion trajectory |15l 
A sketch is provided in Fig. ^ Note that all parameters 
depend on the experimental parameters like materials, 
ion energy and angle of incidence. For amorphous ma- 
terials they can be estimated using the SRIM simulation 
package |2l|. Except where otherwise stated, a — 6.0, 
9 = 50.0°, 4> = 22.0°, and e is chosen to be (2tt) 3 / 2 o-^ 2 , 
while, in contrast to previous studies, we vary systemat- 
ically the values of /i and a mimicing different material 
combinations. 

Furthermore, we extend the sputtering model by con- 
sidering lateral erosion of surface atoms; thus including 
non-linearity in a manner analogous to the Kardar-Parisi- 
Zhang (KPZ) non-linearity that arises from lateral at- 
tachment of adatoms in surface growth [2*2^ . This has 
the effect of relaxing the surface [23j (expecially as local 
surface slopes increase, since h(i) — ► h(i) — rii where n% is 
the number of surface atoms eroded in column i at that 
instant, due to their being closer to the stopping point of 
ion than topmost atom in column (see Fig. 

Surface migration is simulated [161 |24| as a nearest 
neighbor Monte Carlo hopping process, with a site i and 
nearest neighbor site j chosen randomly; and a hop i — ► j 
allowed with probability 

^ocexp(-JL), (2) 

where E = E$ + nnEpj n + Eg b is an energy barrier to 
hopping, comprising of a substrate term Eg = 0.75eV, a 
nearest neighbor term E^n = 0.18eV and a step barrier 
term Esb', nn is the number of in-plane nearest neigh- 
bors of the diffusing atom; T is the substrate tempera- 
ture; and ks is the Boltzmann's constant. Esb = 0.15eV 
if there is no vacant next-nearest neighbors in the plane 
below the hopping atom (at site i), and there is at least 
a vacant next-nearest neighbor in the plane below the 
hopped atom (at site j); otherwise Esb = 0. Thus this 
model discourages diffusion to a down-step edge but once 
the approach is made, it does not discourage a hop down 
the step (i.e since Esb = then; see Fig. 1 of (24|). 
A local heating of the surface occurs right after ion im- 
pact, followed by rapid cooling; hence surface diffusion 
is greatly enhanced, due to higher effective temperatures 
arising from the ion impacts. Consequently, we have used 
a higher effective temperature fc^T = O.leV in our sim- 
ulations, which was estimated in our previous work [H|, 
based on a calculation of the spatio-temporal develop- 
ment of the t emp erature, arising from the local heating 
H| (see also Hf). 

III. RESULTS AND DISCUSSION 



where du = z + a is the distance of the surface atom, from 
the final stopping point of the ion, along the ion trajec- 



The profiles shown in the figures are from sizes L 2 = 
128 , and the bar, on the profiles, indicate the ion beam 
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direction. Distances are in units of a lattice constant, and 
time in ions/atom; except where stated otherwise. By 
varying the collision cascade parameters a and n, we ex- 
plore the full topographic features of this sputter-erosion 
model, for typical ion energies on the order of IKeV; 
keeping a and 9 constant. 

We start with a sketch (Fig. |5J) of the six topographic 
regions that were found, coresponding to different com- 
binations of a and /i, for t — 3 ions/surface atom at 
which almost all the surface topographic features arc dis- 
tinct; the corresponding profiles are shown in Fig|3J Note 
that the boundaries shown in this sketch do not repre- 
sent abrupt transitions from one topography to another, 
as we shall see below (in Figs. ^2 II 01) . Rather we ob- 
serve often a smooth crossover from one behavior to the 
other. Note that the system is still not in (dynamic) equi- 
librium 0], hence the surface morphology still evolves 
with time and at later time a diagram of the type of Fig. 
13 will look different. We have picked a typical time, cor- 
responding to timescales easily accessible in experiments, 
which exhibit a rich behavior as a function of the strag- 
gling parameters fi and a. Also, although the sketch is 
specifically at 9 = 50°, similar "phases" also occur at 
other values of 9, with slight deviations at the bound- 
aries. Examples of a few experimental parameters (ion 
energy, type of material, type of ion etc.) which lead 
to the specific parameters used here we obtained from 
a SRIM simulation [2l|], are; region V: 1.5 - 1.7 KeV 
neon(Ne)-ion, or 2.25 — 2.5KeV argon(Ar)-ion, sputter- 
ing of copper(Cu), or silver(Ag); 1.2 — lAKeV Ne-ion, 
or 1.7— 2.0KeV Ar-ion, sputtering of germanium(Ge), 
or gallium arsenide(GaAs); 170— 200eV helium(He)-ion 
sputtering of graphite(C); region IV: 650 — 800eV Ne-ion 
sputtering of silicon(Si); region III: 800eV — l.lKeV Ar- 
ion sputtering of silicon(Si); 550 — 700eV Ne-ion sput- 
tering of C; with fluences ~ 10 14 — ~ 10 15 ions/cm 2 . 
Note that for most materials and parameter combina- 
tions a < jit, hence the region VI might be difficult ot 
access. Also, SRIM simulations reveal that very large a 
and i.e, beyond the values considered here, are im- 
practical, since they can only occur for a higher a. But 
the value of a is itself restricted by the range of ion ener- 
gies that lead to ripple formation. We are not aware of a 
experimental study of the sputtering behavior were the 
parameters are varied systematically in the fi-a plane. 
Hence, when using these parameters in experiments, one 
might be able to observe the topographies of Fig. |3J 
where some of them have not been observed so far. 

A brief description of each topographic region in 
Fig. |3 including the behavior at later times t > 3 
(atoms/surface atom), is as follows: Region I: rough sur- 
face [see Fig. EJa)] which, as time increases, evolves to a 
hole topography. The sizes of the holes grow and finally 
coalescent to a ripple topography at long times (Fig. 0}. 

Region II: holes are already prominent in this region 
[see Fig. E^b)]; Here the number of holes in this region 
increases with time, and again ripples are formed at long 
times, but at an earlier time than as in region I (not 
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FIG. 2: Different topographic regions (6 = 50°, a = 6). 
Region I: rough surface; II: holes; III: clear ripples oriented 
perpendicular to ion beam direction; IV: short ripples (result- 
ing from increased mu); V: dots; VI: non-oriented structures. 
The arrows indicate the directions referred to in Fig. |H] 




FIG. 3: Profiles for parameters chosen from each topographic 
region in Fig. H 9 = 50°, a = 6, t = 3.0. (I) a = 1, fj, = 0.5; 
(II) a = 1, n = 1.5; (III) a = 3, fi = 1.5; (IV) a = 4, fi = 2.5; 
(V) cr = 5, n = 5; (VI) a = 0.5, (M = 5 



shown as separate figure). 

The number of holes decreases when increasing the 
sputterin depth a [Fig. E|(a)-(c)]. The number of holes 
also decreases with decreasing 9 [Fig. 03 (d)-(f)], while 
ripples are formed already at this early time in this re- 
gion if 9 is increased beyond 9 r (t = 3) ~ 60°. 

Region III: the ripple phase 0,0. Having observed 
in regions I and II that holes evolve into ripples with time, 
we studied this region from the very earliest times (t=0- 
3) but found only very tiny holes, i.e not as pronounced as 
in region II, in the course of ripple formation, see Fig. |HJ). 
9 r ~ 30° in this region. This means, comparing regions 
I, II and three, there seem to be two different processes 
of ripple formation. Ripples can be formed quickly by 
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FIG. 4: Surface profiles of region I, a = 1, // = 1. From 
top-bottom, left-right, * = 3, 20, 40, and 90. 




FIG. 5: Surface profiles of region II, t=3 ions/atom [a — 1, 
H = 1.5). Top row; 9 = 50°, a = 5 (a), a = 6 (b), and a = 7 
(c). Bottom row; a = 6, = 40° (d), 45° (e), and 60° (f). 



evolving directly from a slightly rough surface, or they 
can be formed slowly via the creation of holes, which 
coalescent to ripples on longer time scales. Note that 
in regions I, II, the resulting ripple wavelength is smaller 
than the size of the holes generated on smaller time, while 
in region III the ripple wavelength is larger than the tiny 
holes. 
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FIG. 6: Surface profiles of region III, at very early times. 
From left-right, t=0.1, 0.2 and 0.5 ions/atom (cr = 3, fj. = 1.5, 
9 = 50°, a = 6). 



We also studied here the case of smaller angles, e.g. 



9 = 30°, see Fig. [7] Hence, for lower angles of inci- 
dence, ripple formation is shifted to later times, but the 
ripples tendency again increases when increasing the lon- 
gitudinal struggle a. This indicates that there might be 
a lower critical angle 9 Cl below which no ripple formation 
happens even at long times. Such an effect has been ob- 
served in experiments [9j of sputtering Xe + on Si, where 
below 9c = 40° no ripples have been found at finite but 
long ion fluences. For even smaller angles like 9 = 20°, we 
indeed do not observe ripple formation within the time 
scales (i.e. fluences) we can reach in our simulations. 
Note that a general statement about the existence of such 
a critical angle 6 c (fi,a) would require simulations up to 
very large times for all parameters studied here, which is 
beyond the numerical capacities. Hence, we remain with 
the statement that our numerical results indicate that 
such critical angles indeed exist, without the possibility 
to determine them precisely. 




FIG. 7: 9 = 30°, a = 6, fj, = 1.5; a = 2 (a); a = 3 (b); and 
a = 5 (c). 



Region IV: consists of a mixture of dots and short rip- 
ples, which eventually give way to the dot "phase" (re- 
gion V), as o is increased. Hence, this region seems to 
"interpolate" between regions III and V. 

Region V: consists of dots. These dots are formed on 
some ripple-like structures oriented perpendicular to the 
ion beam direction, as discussed below within more de- 
tail. Noting that our model is a solid on solid model on 
a square lattice, the dots are not unsimilar to the QDs 
predicted by theory pH lllj and observed in experiments 

ED 

Region VI: consists of non-oriented structures exhibit- 
ing a typical lengthscale, but only a slight orientational 
preference parallel to the ion beam. This region, as men- 
tioned above, is probably difficult to access in experi- 
ments. 

A closer look at the dot profiles of region V in Fig. 
reveals the presence of some underlying large-scale 
structures. We now discuss this region and its adja- 
cent regions III and V with more detail. The under- 
lying structure is clearly seen in a 2D structure factor, 
5(k) = (h(k)h(~ k)), where h(k), the fourier transform 
of the height profile h(r), with mean (h) — X) r M r )/^ 2 > 
at time t, is given by 



h(k) 



(h)] exp(ik • r) 



(3) 
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FIG. 8: 2D structure factor (L=128). Top, along ^ = 4 in 
Fig. El left-right: a = 0.5, 1.0, 1.5, and 3.0. Bottom, along 
(7 = 3, left-right: fi = 0.5, 1.0, 1.5, and 3.0 
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FIG. 9: Separation of the 2D S(k) (first image top and bot- 
tom) into regions (in black from second image), for the pur- 
pose of our analysis. The letters I, r, b, denote line, ring and 
background, respectively; and the symbols, _L and ||, denote 
the perpendicular and parallel directions, respectively; as re- 
ferred to in the text. The white regions are excluded from the 
analysis. 



The 2D structure factor shown in Fig. |H1 has been ob- 
tained from an average of 600 independent runs, for pa- 
rameter along [i = 4 (top row), and along a — 3.0 (bot- 
tom row). In these diagrams the structure factor for 
the fc-vektor (k x ,k y ) — (0,0) is displayed in the cen- 
ter and the corners represent the values for (k Xl k y ) — 
27r(±l/8, ±1/8). The case of \x = 4.0, when moving from 
small values of a to larger values (i.e. left-right in Fig. EJ) 
is shown in the top row of Fig. |H1 For small values of /i, 
we are in region IV, where we see a typical wavelength, 
but almost no orientation. This translates to a ring visi- 
ble in the S(k) plot. Note that there is a slight preference 
for an orientation parallel to the ion beam, being visible 
via two peaks in 5(k) at wave vectors perpendicular to 
the ion beam. When increasing a, one moves into region 
V. Here a line perpendicular to the ion beam emerges in 
the structure factor rather abruptly around fi = 1. This 
line represents the underlying structures parallel to the 



ion beam, being visible in in Fig. [3 Note that the "dots" 
emerge on top of these structures; in the 2D structure 
factor their signal is too weak to be visible. 

On the other hand, along a — 3.0, we initially see (bot- 
tom line from bottom-left Fig. [2J an orientation, spread 
around ion beam direction (for instance, at // = 1.0), with 
a more restricted range of k which is typical of the thin 
wobbly ripples. As (i increases, we move to the region 

V as discussed above, but we do not observe an abrupt 
change, because for a large range of values of /i, ripples, 
dots and the underlying parallel structures coexists (re- 
gion III). 

To study the crossover or transition from one region 
to the other in a quantitative way, it might be more in- 
structive to look at order parameters which are numbers 
rahther than the full 2D structure factor. We first define 
a quantity Q = S m (k L )/S m (k||) l where 5 m (k| K _L)) is the 
maximum, for k parallel (perpendicular) to ion beam di- 
rection. This quantity detects the change of orientational 
order. 

When moving along fj, = 4 from region VI into region 

V we expect a rather abrupt change of the behavior from 
the visual inspection of the 2D structure factor (top row 
of |SJ. The behavior of Q when changing a for [i = 4 
is displayed in the inset of Fig. ^| For small values of 
<t, i.e. in region VI, there is only a slight preference of 
structures parallel to the ion beam (corresponding to kj_) 
which leads to small values Q ss 2. When going beyond 
a = 1, i.e. when moving into region V, the domination of 
the longitudinal structures increases, leading to a growth 
ofQ. 

We furthermore want to go beyond studying the height 
of peaks of S(k) by examining the presence or absence 
of a typical length scale in the system. This length scale 
is visible in the 2D structure factor via bright spots (ori- 
ented) or a bright ring. For this purpose, we also look at 
the relative weights 



Qx/y 



J. S(k)da' 



(4) 



of certain areas in the 2D structure factor, where the 
integration da. extends over areas A x /A y of k-points. The 
indices x,y = l,r,b, refer to k values on the line, ring 
and background (outside line and ring) respectively, as 
shown in the upper part of Fig. |5J The white regions, 
in this figure, are excluded from our analysis, because 
they contain a superposition of ks on both the line, and 
the ring. A plot for q r / b is shown in Fig. 1101 other 
combinations (not shown) yield results with information 
not going beyond the study of Q. For q r /f, we observe 
an abrupt change around a = 1 with the ring becoming 
almost indistinguishable {q r /b ~ 1) from the background 
for a > 1 — 1.5, signaling the disappearance of a typical 
length scale around a = 1. Note that q r /jj stays almost 
constant for the full region V, in contrast to the bchvior 
of Q. 

No we turn to the behavior along the a = 3 line, i.e. 
when moving from region I, into region III, then IV and 
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1 2 (j 3 4 5 

FIG. 10: Plot ol the quantities q r /t and Q (inset) as a func- 
tion of a, along fj, = 4.0, as described in the text. 



finally V. The behavior of Q as function of /j, (bottom 
row of Fig. [HJ is displayed in the inset of Fig. ^] The 
crossover around /i ~ 2 from the ripple region (II) where 
k\\ wave vectors dominate {Q < 1) to the region V where 
/c^-structures dominate is clearly visible. 

The same result is obtained, when again studying not 
only peaks, but integrated structure factors over certain 
selected areas. Hence, we also study fc-vectors parallel 
and perpendicular to the ion beam (x,y =|| , _L), as shown 
in the lower part of Fig. [5] Also q±/\\ exhibits a strong 
growth for fi > 2, i.e. when moving into region III. 

Finally, a study of q r (not shown) again confirms the 
the loss of a certain length scale when moving into region 
V. This is again visible via the disappearance of weight 
on the ring in the 2D structure factor with respect to 
the background [see again Fig. [8] . In any case, we have 




1 2 3 4 5 



FIG. 11: Plot of the quantities q±/\\ and Q (inset) as a func- 
tion of fi, along a = 3.0, as described in the text. 

mainly looked at snapshots at time t = 3, correspond- 
ing to typical experimental fluences of 10 15 ions/cm 2 (as 
already stated above). Hence, the location of the exact 
boundaries between the different regions, is not a major 
interest here. We only want to demonstrate, that indeed 
order-parameter-like functions can be defined. The full 
information is anyway contained in the 2D structure fac- 
tors. 



In the asymptotic limit, the morphology of sput- 
tered surfaces may be described by the noisy Kuramoto 
Sivashinsky equation 0] 

8th = v\ l d 2 h+v ± d 2 y h+^(d x h) 2 + ^(d y h) 2 -Kd i h+T 1 , 

(5) 

where v\\{y±), Ay(Aj_), is the effective surface tension co- 
efficient, and nonlinear coefficient, respectively, along a 
direction parallel (perpendicular) to the projection of the 
ion beam direction onto the surface plane (here, the x- 
axis is parallel to the projection); K is the surface dif- 
fusion constant; and rj is an uncorrelated white noise 
with zero mean. These coefficients are explicitly given, 
in terms of the sputtering parameters, in ref. Note 
that erosion of surface material tends to maximize the 
exposed area, so-called negative surface tension, hence, 
for the sputtering phenomena, this instability constrains 
the coefficients and v± to be negative. 

At early times, the local slopes are small enough in 
most of the regions for us to ignore the nonlinearities; we 
are therefore left with a noisy Bradley-Harper equation 
[jj. A plot of the coefficients v for our parameter range, 
along a = 3, is shown in Fig. 1121 |mi | and \v±_ | are nonzero 
always which implies the presence of two lengthscales. 
Hence, according to the linear continuum theory, ripples 
parallel and perpendicular to the ion beam direction are 
always present in the system; the one observed is the one 
for which \v\ is highest; i.e, with the highest amplitude- 
growth rate R = -{v\\k 2 + v±k 2 L ) - K(kj + k\) 2 0. 

We now try to understand roughly the behvior along 
the above discussed a — 3 and fi = 4 lines within the 
linear theory. For a full understanding one would have 
to consider also the nonlinear terms, where the full de- 
pendance on the parameters is yet not available to the 
authors. In Fig. ^|the values of i/y and u± are shown as 
function of fi (main plot) and a (inset) along these lines. 
For the case fx = 4 (see inset), we observe that v\\ > 
for all values of er, hence the preferential orientation is 
always parallel to the ion beam, as observed. To under- 
stand the crossover from region VI to region V, one has 
probably to consider the nonliear terms. 

For the a — 3 line, we observe \v±\ w \vn \ for /i < 1, 
which is compatible with the behavior in region I, where 
no preferential orientation is observed. For 1 < // < 2.5, 
1^11 1 > I^J-I which implies the dominance of ripples with 
fey (i.e, with wavelength A = 2n^/2K/\is^ |), as we have 
seen in region III and partly in the crossover region IV. 
In region IV ripples are still present, but dominate less. 
For fi above this range, the structures with k± dominate. 
However, the positive value of un for higher fj, is contrary 
to our results, since it implies a preferential smoothening 
along the parallel direction. We have not observed such 
smoothening, here probably nonlinear effects are more 
important. If we consider the maximum of \ and \v±\ 
for the region where v\\ < and/or v p erp < 0, then we 
see that this maximum is obtained in region V for large 
fi. Indeed we observe that region V is rougher than e.g. 
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FIG. 12: The surface tension coefficients of Eq. |SJ along 
a = 3.0, for our simulation parameters. 6 — 50°, a — 6.0 




FIG. 14: Surface width W versus time, a = 5, /i = 5 (region 
V), 8 = 50°. In the inset, time is in multiples of 10 2 ; there are 
two scaling regions t = 40 — 140 and t = 1000 — 1800, which 
yields /3 = 0.665 and 0.137 respectively. 



region III, see below. 



3= 4 







— P=0-327 




A 0.625 


1 1 1 1 1 1 1 1 1 i .. 


/ ^ 0.5 






□ 


' 0.375 
0,3 




n □□□□^ifflmmiiiftir 
i i 




20 30 40 50 t 80100 150 



t 

FIG. 13: Surface width W versus time, a = 3, fj, = 1.5 (region 
III), 9 = 50°. In the inset is a plot of the consecutive /3s for 
t = 20 — 200, where the lower and upper cutoffs are seen to 
occur at t = 35 and 150 respectively. 

But, before we proceed, we should recall that regions 
I and II eventually evolve to a similar ripple topogra- 
phy as in region III; region IV is an intermediate stage 
and VI is hardly experimentally accessible. Therefore, 
we are interested in two different scaling regions, i.e, aris- 
ing from ripple-yielding, as well as dot yielding parame- 
ters. We use the finite-size scaling behavior of the surface 
width/roughness, 

W(L,t) = L<*f(±), (6) 

where a and (5 are the roughness exponent and growth 
exponent, respectively; and W — ^ r {h{r) — (h)) 2 /L 2 , 
i.e, W is the root-mean-squared fluctuations in the sur- 
face height. I@ defines two scaling regions, separated at 
the saturation time t s (L) oc L z . z — a/ [3 is the dynamic 



exponent. In the first region (t << L z ) W oc t@ , indepen- 
dent of L; and in the second region (t >> L z ), W cx L a , 
constant in time [23|. A scaling argument in reveals 
that t s increases with the diffusion coefficient, and with 
decreasing strength of non-linearities 28] . 

In the following roughness analysis, we rescale the time 
unit such that L 2 particles are eroded in unit time, which 
is analogous to the measure of time in simulations of 
epitaxial growth on vicinal surfaces. 

We obtained our scaling exponents from an average of 
600 independent runs. Fig. ED is a plot of the surface 
width W(t) for L — 128, we determine the actual scaling 
region of (3 from the consecutive slopes [2^ of W(t) for 
t < L z ; shown in the inset of Figs. 1 1 31 and PHI where a fit 
to this region gives (3 = 0.327 ± 0.001 and j3 = 0.665 ± 
0.003 (from the scaling for times around t = 100) for the 
topographic regions (see Fig. 0) III and V respectively. 

Our exponents are quite different from the scaling of 
the KPZ equation ^22j in (2+1) dimensions [13, |29|| . 
which shows that for these parameter choices, either one 
or both of the negative surface tension coefficients, in the 
noisy KS equation 14], do not renormalize to positive 
values. That is, if indeed, this atomistic model which 
is known to accurately describe the morphology of sput- 
tered surfaces at early and long times |l5Lll6Ll30j ]. can be 
represented by the noisy KS equation, in the asymptotic 
limit, as expected for the sputtering phenomenon. 



IV. CONCLUSIONS 

We have extended the study of surface topogra- 
phy induced by ion bombardment beyond the circular- 
symmetric collision cascade case a = fj, = 1; using a 
discrete Monte Carlo model. We found an upper bound 
fi ~ 2 for clear ripple formation. We observed a crossover 
from hole topography to ripple topography, for same col- 
lision cascade parameters. Above the upper bound we 
found a crossover to a dot topography; our results in- 



8 



dicate the non-trivial influence of the value of [i for the 
surface topography. Moreover, we found the possibility 
of dot formation for \i > 2 and a > 1.7, without sample 
rotation. 

Also, we found the possibility of different incidence 6 r , 
around which ripple formation is possible, as in the ex- 
periment; our results indicating that high value of the 
longitudinal straggling result in lower value of 8 r . Finally, 
we found that the collision cascade parameters affect the 
growth exponent. 

Our results indicate that using new target/projectil 
combinations different in experiments might lead to sur- 
face topographies yet not observed. Here a systematic 
study covering the full parameter space (at least in the 
range (<r < n) would be very interesting, which has not 
been done so far to our knowledge. Furthermore, it would 
be nice to investigate the small-time behavior, i.e. the rip- 
ple formation process, experimentally. In this way one 



could verify whether the two different creation mech- 
anism we have observed, hole-coalescence and creation 
from a rough surface, can be observed in experimental 
systems. 
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